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I.  Introduction 

The  research  supported  by  this  grant  centered  on  the  design  of  theoretical  and 
computational  approaches  to  predict  promising  High  Energy  Density  Matter  (HEDM) 
and  to  guide  the  efficient  synthesis  of  these  materials.  The  first  objective  was  to  develop 
methodology  for  the  simulation  of  photochemical  (i.e.,  electronically  excited)  reactions, 
which  could  play  an  important  role  in  the  design  and  synthesis  of  HEDM.  The  second 
objective  was  to  develop  methodology  for  quantum  dynamical  calculations  of  polyatomic 
reactions  to  aid  in  the  prediction  and  synthesis  of  new  energetic  materials.  The  third 
objective  was  to  investigate  the  solvent  effects  for  fundamental  organic  reactions  to  aid  in 
the  efficient  synthesis  of  HEDM.  The  fourth  objective  was  to  develop  methodology  for 
the  calculation  of  hydrogen  vibrational  wavefiinctions  to  enable  the  simulation  of 
hydrogen  transfer  reactions  required  for  the  synthesis  of  polyhedral  oligomeric 
silsesquioxanes  (POSS),  which  are  used  in  coatings  that  are  resistant  to  extreme 
conditions  such  as  heat  and  abrasion. 

II.  Methods  and  Results 

A.  Simulation  of  electronically  excited  reactions  . 

We  have  developed  new  methodology  for  the  simulation  of  photoexcited 
reactions.  In  these  reactions,  absorbed  light  promotes  the  system  to  a-highly  reactive 
electronically  excited  state,  typically  leading  to  products  of  higher  energy  than  the 
reactants.  Thus,  photochemical  reactions  could  play  an  important  role  in  the  design  and 
synthesis  of  HEDM.  One  example  is  the  photochemical  ring  closure  of  a  polycyclic 
dialkene  with  two  double  bonds  that  react  to  produce  a  highly  strained  cage  molecule 
such  as  cubane  (C6H8),  which  is  a  promising  HEDM  candidate.  The  ability  to 
theoretically  study  the  dynamics  of  photochemical  reactions  under  a  range  of 
experimental  conditions  will  help  guide  the  efficient  synthesis  of  promising  HEDM 
candidates. 

The  simulation  of  photoexcited  reactions  is  challenging  since  multiple  coupled 
electronic  states  are  involved.  Thus,  nonadiabatic  quantum  mechanical  effects  must  be 
incorporated  into  the  simulations.  Fully  quantum  mechanical  simulations  are 
impractical  for  systems  containing  more  than  a  few  atoms.  As  a  result,  mixed 
quantum/classical  methods  have  been  developed  for  the  simulation  of  larger  systems. 

For  photoexcited  reactions,  typically  the  quantum  subsystem  consists  of  the  electrons, 
and  the  classical  subsystem  consists  of  the  nuclei.  One  particularly  successful  approach 
for  simulating  photoexcited  reactions  is  traj  ectory  surface  hopping.  Our  research  has 
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focused  on  the  Molecular  Dynamics  with  Quantum  Transitions  (MDQT)  surface  hopping 
method.1'3  In  MDQT  the  classical  subsystem  is  approximated  as  an  ensemble  of 
independent  trajectories,  and  each  trajectory  moves  classically  on  a  single  potential 
energy  surface  with  the  possibility  of  instantaneous  transitions  among  the  surfaces.  The 
quantum  amplitudes  for  all  surfaces  are  propagated  coherently  along  each  independent 
trajectory,  and  the  probability  of  a  transition  depends  on  the  rate  of  change  of  the 
quantum  probabilities  determined  from  the  quantum  amplitudes.  The  number  of 
transitions  is  minimized  by  specifying  that  the  flux  of  trajectories  switching  from  one 
state  to  another  is  uni-directional  over  a  specified  time  interval.  This  algorithm  is 
designed  to  ensure  that  the  fraction  of  trajectories  on  each  surface  is  equivalent  to  the 
corresponding  average  quantum  probability.  As  has  been  noticed  in  the  literature, 
however,  this  internal  consistency  is  not  always  maintained. 

The  specific  aim  of  this  project  was  to  identify  the  reasons  for  this  discrepancy 
and  to  develop  methods  for  improving  the  internal  consistency  of  MDQT.  We  have 
identified  two  reasons  for  the  internal  inconsistency  in  trajectory  surface  hopping 
methods.  These  reasons  are  the  existence  of  classically  forbidden  transitions  and  the 
divergence  of  the  independent  trajectories.  We  have  developed  a  modified  MDQT 
method  that  improves  the  internal  consistency.4  In  this  method,  the  classically  forbidden 
transitions  are  eliminated  by  setting  the  nonadiabatic  coupling  between  each  pair  of  states 
to  zero  for  the  integration  of  the  time-dependent  Schrodinger  equation  if  the 
corresponding  transition  would  be  classically  forbidden .  This  approach  eliminates  the 
classically  forbidden  transitions  since  the  probability  of  a  quantum  transition  between 
two  states  vanishes  if  the  nonadiabatic  coupling  between  the  two  states  is  zero. 

(Note  that  this  procedure  completely  avoids  the  controversial  issue  of  velocity  rescaling 
after  a  classically  forbidden  transition.)  In  addition  to  this  modification,  the  difficulties 
due  to  divergent  trajectories  are  alleviated  by  removing  the  coherence  of  the  quantum 
amplitudes  when  each  trajectory  leaves  a  nonadiabatic  coupling  region. 

We  have  compared  the  standard  and  modified  MDQT  methods  to  fully  quantum 
calculations  for  a  classic  model  for  ultrafast  electronic  relaxation  (i.e.,  a  two-state  three¬ 
mode  model  of  the  conically  intersecting  Si  and  S2  excited  states  of  pyrazine) 4  While 
the  standard  MDQT  calculations  exhibit  significant  discrepancies  between  the  fraction  of 
trajectories  in  each  state  and  the  corresponding  average  quantum  probability,  the 
modified  MDQT  method  leads  to  remarkable  internal  consistency  for  this  model  system. 
For  this  model,  the  divergence  of  trajectories  is  mainly  responsible  for  internal 
consistency,  but  the  elimination  of  the  classically  forbidden  transitions  also  leads  to  a 
minor  improvement.  This  methodology  will  be  useful  for  the  simulation  of  photoexcited 
reactions  aimed  at  the  synthesis  of  promising  HEDM  candidates. 

B.  Development  of  TDSCF-RPH  methodology 

We  have  developed  new  methodology  that  allows  quantum  dynamical 
calculations  of  general  polyatomic  reactions.  This  methodology  will  aid  in  the  prediction 
and  synthesis  of  new  energetic  materials.  We  developed  the  TDSCF-RPH  (time- 
dependent  self-consistent-field  reaction  path  Hamiltonian)  dynamics  method,  which 
allows  the  calculation  of  the  real-time  quantum  dynamics  of  chemical  reactions  involving 
polyatomic  molecules.5  7  This  approach  is  based  on  a  reaction  path  Hamiltonian  (RPH), 
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which  characterizes  the  reaction  dynamics  in  terms  of  motion  along  a  reaction  path  and 
vibrational  motion  orthogonal  to  this  reaction  path.  In  this  approach  the  reaction  is 
assumed  to  occur  within  a  multidimensional  harmonic  valley  about  the  reaction  path. 

The  main  advantage  of  such  RPHs  is  that  they  require  only  a  relatively  small  number  of 
ab  initio  calculations  of  the  potential  energy  surface. 

Reaction  path  Hamiltonians  have  been  derived  for  different  types  of  reaction 
paths.  Typically  the  reaction  path  is  assumed  to  be  the  minimum  energy  path  (MEP), 
which  is  defined  as  the  path  of  steepest  descent  in  mass-weighted  coordinates  starting 
from  a  transition  state  down  toward  reactants  and  toward  products.  Miller,  Handy,  and 
Adams  derived  a  RPH  that  describes  the  reaction  dynamics  of  a  polyatomic  system  in 
terms  of  motion  along  the  MEP  and  vibrational  motion  orthogonal  to  the  MEP.  This 
MEP  RPH  has  been  discussed  extensively  in  the  literature  and  has  been  used  to  describe 
a  variety  of  chemical  processes.  Although  highly  successful  in  many  areas,  this  MEP 
RPH  does  not  provide  an  appropriate  description  of  the  reaction  dynamics  for  situations 
in  which  the  reaction  path  is  sharply  curved  since,  in  this  case,  the  dynamical  motion 
deviates  far  from  the  MEP.  To  account  for  such  situations.  Miller  and  coworkers  have 
developed  a  diabatic  straight-line  RPH  that  interpolates  linearly  between  equilibrium 
reactant  and  product  geometries  of  the  molecular  system. 

We  have  developed  the  theory  and  numerical  methodology  for  TDSCF  quantum 
dynamics  based  on  the  reaction  path  scheme.  Our  first  TDSCF-RPH  method  was  based 
on  the  MEP  RPH  of  Miller,  Handy,  and  Adams.  We  derived  the  TDSCF  equations  of 
motion  for  cases  in  which  the  MEP  has  zero  or  small  reaction  path  curvature  and  showed 
that  when  the  coupling  between  the  normal  modes  is  negligible  the  dynamics  reduces  to  a 
one-dimensional  numerical  time  propagation.5  We  also  numerically  tested  this  TDSCF- 
RPH  methodology  by  comparison  to  exact  quantum  dynamics  based  on  the  exact 
Hamiltonian  for  a  series  of  simple  model  systems.6  The  success  of  this  methodology  for 
these  simple  model  systems  is  encouraging.  ,  ^ 

Our  second  TDSCF-RPH  method  was  based  on  diabatic  RPHs.7  The  diabatic 
representation  of  the  RPH  involves  a  transformation  of  variables  that  shifts  the  coupling 
between  the  vibrational  modes  from  the  kinetic  energy  to  the  potential  energy.  The 
methodology  we  developed  is  applicable  to  the  MEP  RPH  for  the  case  in  which  the  MEP 
has  negligible  curvature  and  to  the  straight-line  RPH  for  the  case  in  which  the  MEP  has 
large  curvature  but  the  system  follows  a  straight-line  path  instead  of  the  MEP.  We 
derived  the  equations  of  motion  for  diabatic  RPHs  and  showed  that  in  the  diabatic 
representation  the  TDSCF  dynamics  reduces  to  a  one-dimensional  numerical  time 
propagation,  even  for  the  case  of  large  coupling  between  the  normal  modes.  These 
numerically  efficient  approaches  for  simulating  polyatomic  reactions  will  aid  in  the 
prediction  and  synthesis  of  HEDM. 

C.  Solvent  effects  on  organic  reactions 

We  have  investigated  the  solvent  effects  for  fundamental  organic  reactions.  The 
two  reactions  we  have  studied  are  the  Claisen  rearrangement  of  allyl  vinyl  ether  and  the 
Diels  Alder  reaction  of  cyclopentadiene  with  methyl  vinyl  ketone.8  Both  of  these 
reactions  are  useful  in  synthesis,  so  elucidation  of  solvent  effects  could  aid  in  the  efficient 
synthesis  of  HEDM.  The  methodology  developed  is  also  applicable  to  other  reactions. 
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The  reaction  rates  of  both  Claisen  rearrangements  and  Diels  Alder  reactions  have 
been  found  to  be  accelerated  in  polar  solvents.  Experimental  data  suggest  that  the 
Claisen  rearrangement  of  allyl  vinyl  ether  is  accelerated  by  as  much  as  a  factor  of  1000  in 
water  relative  to  the  gas  phase.  Moreover,  experiments  also  indicate  that  the  Diels  Alder 
reaction  of  cyclopentadiene  with  methyl  vinyl  ketone  is  accelerated  by  a  factor  of  730  in 
water  relative  to  isooctane. 

Both  the  Claisen  rearrangement  of  allyl  vinyl  ether  and  the  Diels  Alder  reaction  of 
cyclopentadiene  with  methyl  vinyl  ketone  have  been  studied  with  a  wide  range  of 
computational  methods.  One  approach  has  been  to  perform  electronic  structure 
calculations  on  the  reacting  system  in  the  presence  of  a  dielectric  continuum  solvent. 
Typically  these  calculations  are  problematic  due  to  the  lack  of  consideration  of  explicit 
hydrogen  bonding  interactions.  An  alternative  approach  has  been  to  perform  Monte 
Carlo  or  molecular  dynamics  simulations  for  the  reacting  system  in  a  periodic  box 
containing  hundreds  of  explicit  solvent  molecules.  These  calculations  have  been 
successful  in  reproducing  experimentally  measured  free  energies  of  solvation.  These 
calculations  are  extremely  computationally  intensive,  however,  and  most  studies  neglect 
the  dynamical  solvent  effects.  (One  notable  exception  is  the  study  by  Voth  and 
coworkers,  which  includes  a  reactive  flux  calculation  of  the  dynamical  transmission 
coefficient.) 

We  have  investigated  the  solvent  effects  for  these  two  reactions  within  the 
framework  of  the  reaction  path  Hamiltonian  (RPH)  developed  by  Miller,  Handy,  and 
Adams.8  This  RPH  is  based  on  the  minimum  energy  path  (MEP)  described  above.  The 
initial  step  of  this  project  was  the  generation  of  the  MEPs  at  the  DFT/B3LYP/6-31G* 
level  for  these  reactions  in  the  gas  phase  and  in  the  presence  of  two  explicit  water 
molecules.  We  found  that  the  presence  of  the  two  water  molecules  lowered  the  activation 
energy  barrier  due  to  stronger  hydrogen  bonding  between  the  water  molecules  and  the 
solute  oxygen  in  the  transition  state  than  in  the  reactant.  The  magnitude  of  the  calculated 
free  energy  of  solvation  agreed  qualitatively  with  experimental  estimates  and  other 
calculations.  We  also  analyzed  the  geometries  and  partial  charges  along  the  MEPs  to 
determine  the  structural  and  electrostatic  roles  of  the  water  molecules.  In  both  cases  we 
found  that  the  presence  of  the  water  molecules  resulted  in  a  looser,  more  dissociative 
transition  state.  We  also  found  that  the  transition  state  was  more  polar  in  the  presence  of 
the  water  molecules. 

In  addition,  we  analyzed  the  frequencies,  couplings,  and  curvatures  along  the 
MEPs  in  terms  of  the  RPH  for  the  two  reactions  in  the  gas  phase  and  in  the  presence  of 
two  water  molecules.8  The  Claisen  rearrangement  exhibits  a  curvature  peak  on  the 
reactant  and  the  product  side,  but  the  Diels  Alder  reaction  exhibits  only  a  peak  on  the 
product  side.  For  both  reactions,  the  overall  curvature  increases  in  the  presence  of  the 
two  water  molecules  due  to  coupling  of  the  solvent  to  the  reaction  coordinate.  The 
strongly  coupled  modes  in  the  gas  phase  contain  contributions  from  the  stretching  motion 
of  bonds  that  are  broken  and  formed  during  the  reaction.  In  the  presence  of  two  water 
molecules,  the  strongly  coupled  modes  contain  additional  contributions  from  the 
stretching  motion  of  the  hydrogen  bonds  between  the  water  molecules  and  the  solute 
oxygen.  We  found  that  the  solvation  increases  the  curvature  more  for  the  Diels  Alder 
reaction  than  for  the  Claisen  rearrangement,  indicating  that  the  water  molecules  are  more 
strongly  coupled  to  the  reaction  coordinate  for  the  Diels  Alder  reaction. 
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We  have  calculated  the  dynamical  transmission  coefficients  for  both  reactions 
using  reactive  flux  molecular  dynamics  methods  based  on  a  reaction  path  Hamiltonian. 
Due  to  the  locations  of  the  coupling  peaks,  reflections  occur  in  both  the  reactant  and  the 
product  regions  for  the  Claisen  rearrangement  but  only  in  the  product  region  for  the  Diels 
Alder  reaction.  (For  both  reactions,  reflections  may  also  occur  in  the  transition  state 
region.)  Recrossings  of  the  transition  state  in  the  forward  direction  were  observed  for  the 
Diels  Alder  reaction  but  not  for  the  Claisen  rearrangement.  The  transmission  coefficients 
were  nearly  unity  for  the  Claisen  rearrangement  both  in  the  gas  phase  and  in  the  presence 
of  two  water  molecules.  In  contrast,  the  transmission  coefficients  were  slightly  smaller 
for  the  Diels  Alder  reaction:  k=0.95  in  the  gas  phase  and  k=0.89  in  the  presence  of  two 
water  molecules.  The  smaller  transmission  coefficients  for  the  Diels  Alder  reactions  are 
due  mainly  to  the  flatter  potential  energy  along  the  reaction  coordinate  near  the  transition 
state  and  the  closer  proximity  of  the  coupling  peak  to  the  transition  state.  These 
differences  lead  to  lower  momentum  along  the  reaction  coordinate  in  the  strong  coupling 
region  and  thus  a  higher  probability  of  reflections.  The  decrease  in  the  transmission 
coefficient  for  the  Diels  Alder  reaction  in  the  presence  of  two  water  molecules  is  due  to  a 
combination  of  a  flatter  potential  energy  along  the  reaction  coordinate  near  the  transition 
state  and  a  larger  curvature  peak  on  the  product  side. 

Our  investigation  of  the  solvent  effects  for  these  two  fundamental  organic 
reactions  elucidated  the  electrostatic,  structural,  and  dynamical  roles  of  water  molecules, 
as  well  as  possible  reaction  mechanisms.  This  approach  is  less  computationally  intensive 
than  full-scale  molecular  dynamics  simulations  but  still  provides  useful  qualitative 
information  about  dynamical  solvent  effects.  Moreover,  this  general  approach  can  be 
used  to  study  the  effects  of  different  solvents.  This  type  of  comparative  study  will  be 
useful  for  the  efficient  synthesis  of  HEDM  candidates. 

D.  Calculation  of  hydrogen  vibrational  wavefunctions  ** 

We  have  developed  new  methodology  for  the  calculation  of  hydrogen  vibrational 
wavefunctions  for  systems  in  which  only  one  or  a  few  specified  hydrogen  atoms  are 
treated  quantum  mechanically,  and  the  remaining  nuclei  remain  fixed.9'  This 
methodology  will  enable  the  simulation  of  hydrogen  transfer  reactions  required  for  the 
synthesis  of  polyhedral  oligomeric  silsesquioxanes  (POSS).  POSS  is  important  to  the 
HEDM  project  for  use  in  coatings  that  are  resistant  to  extreme  conditions  such  as  heat 
and  abrasion. 

The  methodology  we  developed  to  calculate  hydrogen  vibrational  wavefunctions 
is  called  the  FGH-MCSCF  (Fourier  grid  Hamiltonian  multiconfigurational  self- 
consistent-field)  method.9  In  the  FGH-MCSCF  method,  the  potential  energy  surface  for 
the  transferring  hydrogen(s)  is  calculated  on  a  grid  (i.e.,  all  nuclei  except  the  transferring 
hydrogen(s)  are  fixed,  and  the  energy  is  calculated  for  the  hydrogen(s)  positioned  at  each 
point  on  the  grid).  The  hydrogen  vibrational  wavefunctions  are  calculated  by 
numerically  solving  the  time-independent  Schrodinger  equation  for  the  hydrogen  nucleus 
(or  nuclei)  moving  on  this  potential  energy  surface.  The  multidimensional  wavefunctions 
are  calculated  with  an  MCSCF  approach,  where  each  hydrogen  vibrational  wavefunction 
is  expressed  as  a  linear  combination  of  single  configurations  that  are  products  of  one¬ 
dimensional  wavefunctions  represented  directly  on  the  grid.  The  variational  method  is 
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utilized  to  minimize  the  total  energy  with  respect  to  both  the  configuration  coefficients 
and  the  one-dimensional  wavefunctions.  We  have  found  that  only  a  relatively  small 
number  of  configurations  are  required  for  the  calculation  of  accurate  hydrogen 
vibrational  wavefunctions.  Thus,  the  FGH-MCSCF  approach  avoids  the  diagonalization 
of  large  matrices  required  for  a  full  configuration  interaction  approach.  In  order  to 
further  increase  the  efficiency  of  the  calculation  of  nuclear  wavefunctions,  we  have 
developed  a  partial  multidimensional  grid  generation  method.11  This  method 
substantially  decreases  the  number  of  potential  energy  calculations  by  avoiding  this 
calculation  for  grid  points  with  high  potential  energy.  We  have  tested  both  the  FGH- 
MCSCF  and  the  partial  multidimensional  grid  generation  methods  by  application  to 
hydride  transfer  in  model  systems.10,1 1  This  methodology  will  be  useful  for  the 
investigation  of  hydrogen  transfer  reactions  required  for  the  synthesis  of  POSS. 

III.  Summary 

Thus,  we  have  achieved  the  four  research  objectives  listed  in  the  Introduction. 
First,  we  further  developed  the  Molecular  Dynamics  with  Quantum  Transitions  (MDQT) 
surface  hopping  method  for  the  simulation  of  electronically  excited  reactions.  The 
internal  consistency  of  MDQT  was  improved  by  eliminating  classically  forbidden 
transitions  and  by  damping  the  coherence  of  the  quantum  amplitudes.  Second,  we 
developed  the  TDSCF-RPH  (time-dependent  self-consistent-field  reaction  path 
Hamiltonian)  dynamics  method  for  the  calculation  of  the  real-time  quantum  dynamics  of 
chemical  reactions  involving  polyatomic  reactions.  In  addition  to  the  development  of  the 
theory  and  computational  algorithm  for  TDSCF-RPH,  numerical  tests  were  performed  on 
model  systems.  Third,  we  investigated  the  solvent  effects  for  a  Claisen  rearrangement 
and  a  Diels  Alder  reaction.  Our  approach  combined  electronic  structure  calculations  with 
reactive  flux  molecular  dynamics  calculations  based  on  a  reaction  path  Hamiltonian. 
Analysis  of  the  dynamical  trajectories  provided  insight  into  the  dynamical  role  of  the 
water  molecules  and  elucidated  possible  reaction  mechanisms.  Fourth,  we  developed  the 
FGH-MCSCF  (Fourier  grid  Hamiltonian  multiconfigurational  self-consistent-field)  and 
partial  multidimensional  grid  generation  methods  for  the  efficient  calculation  of 
multidimensional  hydrogen  vibrational  wavefunctions.  These  two  methods  were  applied 
to  hydride  transfer  in  a  model  system.  All  of  these  new  theoretical  and  computational 
approaches  will  aid  in  the  efforts  of  the  Air  Force  to  predict  promising  energetic  materials 
and  to  guide  the  efficient  synthesis  of  these  materials. 
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